Mixed-Effects Models
Disturbance through human hunting activity can significantly impact prey species through both consumptive and nonconsumptive effects. The nonconsumptive effects of rabbit hunting on Northern Bobwhite (Colinus virginianus; hereafter, bobwhite) could cause an increased perceived risk of predation by bobwhite during rabbit hunting events may elicit anti-predator responses, such as reduced movement away from the safety of cover.
bobwhite <- read.csv('bobwhite3.csv')
bobwhite$ID <- as.factor(bobwhite$ID) #make ID a factor
p1<- ggplot(bobwhite, aes(x=HuntDay, y=HW_Dist, group=ID, color=ID, shape=ID)) +
geom_point(size=4, alpha=0.6, position = position_dodge2(width=.33, preserve = "total")) +
scale_y_continuous() +
#geom_line() +
geom_smooth(method = "lm", se = FALSE) +
labs(title="Risk Behavior in Bobwhite During Hunting Season", x= "Hunting Season Species", y = "Distance from Hardwood Forest Cover (meters)")+
theme_bw()+
scale_color_brewer(palette = "BrBG")
p1## `geom_smooth()` using formula 'y ~ x'
# ID: Unique ID given to each bobwhite covey tracked in chronological order.
# HuntDay: Denotes if it was a “Rabbit” or “Quail” (bobwhite) scheduled hunt day.
# HW_Dist: Distance in meters a bobwhite covey was from hardwood habitat.bobwhite_means <- bobwhite %>%
group_by(HuntDay) %>%
summarise(mean_HW_Dist=mean(HW_Dist),
se_HW_Dist=sd(HW_Dist)/sqrt(n()))
bobwhite_means## # A tibble: 2 × 3
## HuntDay mean_HW_Dist se_HW_Dist
## <chr> <dbl> <dbl>
## 1 Quail 22.3 4.29
## 2 Rabbit 57.0 5.30mixed_bobwhite <- lmer(HW_Dist~(HuntDay*ID)+(1|ID), data = bobwhite)## Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
## eigenvalue close to zero: 5.5e-09anova(mixed_bobwhite) ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## HuntDay 19204.9 19204.9 1 188 11.0107 0.001088 **
## ID 3168.5 633.7 5 188 0.3633 0.873151
## HuntDay:ID 23341.2 4668.2 5 188 2.6764 0.023103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1performance::check_model(mixed_bobwhite)
bob_means <- emmeans(mixed_bobwhite, "HuntDay")## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installed## NOTE: Results may be misleading due to involvement in interactionsbob_means## HuntDay emmean SE df lower.CL upper.CL
## Quail 23.1 25.3 188 -26.92 73.1
## Rabbit 45.1 25.1 188 -4.31 94.5
##
## Results are averaged over the levels of: ID
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95Because there is an unequal number of points per day for each individual covey, we can check the accuracy of this model by using emmeans:
em_bob <- as.data.frame(bob_means)
em_bob## HuntDay emmean SE df lower.CL upper.CL
## 1 Quail 23.08173 25.34776 188 -26.920856 73.08431
## 2 Rabbit 45.11574 25.05572 188 -4.310749 94.54224bobwhite_means## # A tibble: 2 × 3
## HuntDay mean_HW_Dist se_HW_Dist
## <chr> <dbl> <dbl>
## 1 Quail 22.3 4.29
## 2 Rabbit 57.0 5.30em_bob$HuntDay <- factor(em_bob$HuntDay, levels= c("Quail","Rabbit"))
ggplot(em_bob, aes(x=HuntDay, y=emmean)) +
geom_point(size=5, color="#87b8d3") +
geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), width=.2, color="#6699CC") +
geom_point(size=5, data=bobwhite_means, x=bobwhite_means$HuntDay, y=bobwhite_means$mean_HW_Dist, color = "#336699") +
theme(axis.text.x = ggtext::element_markdown(color = "tan4", size = 12)) +
scale_x_discrete(labels = labels) +
theme(plot.caption=element_text(size=9, hjust=0, margin=margin(15,0,0,0)))+
theme_bw()+
labs(title="Comparing raw means to emmeans", caption="Light blue = raw means, Dark blue = adjusted means", x= "Hunting Day (1=Quail, 2=Rabbit)", y = "Average distance from hardwood cover (m)") 
\(~\) \(~\)
Nested Hierarchical Model
Great Tits (Parus major)
Advertisement signaling is usually linked to intersexual selection and intrasexual competition and thus is a key component of a species’ ecology. Using a novel spatial tracking system, the authors tested whether or not the spatial behavior of male and female great tits (Parus major) changes in relation to the response of a territorial male neighbor to an intruder. They tracked the spatial behavior of male and female great tits (N = 13), 1 hr before and 1 hr after simulating territory intrusions. The individual bird is a random effect, and the sex of the bird and measurements taken before and after are fixed effects.
tit <- read.csv("great_tits.csv")
tit$sex <- as.factor(tit$sex)
ggplot(tit, aes(sex, dist_m, colour = as.factor(id), shape=as.factor(measure))) +
geom_jitter(width =0.15, size=5, alpha=0.6)+
ylab ("Response Distance from Playback Location (meters)") +
xlab ("Sex") +
annotate("text", x = 2, y = 76, label = "13 birds") +
annotate("text", x = 2, y = 72, label = "2 measurements per bird") +
annotate("text", x = 2, y = 67, label = "26 total measurements", size=4, color="blue")+
labs(title="Sex-specific Responses to Territorial Intrusions", caption="1 = female
2 = male")+
scale_shape_discrete(name= "Measurements
Before / After")+
theme_bw()
In this study, 13 individuals were measured 2 times, 1 hour before playback and 1 hour after. The individuals are separated by sex to look at the different distances between females and males. Distance on the Y axis is measured in meters, with distance being how many meters away a bird was detected before and after a playback recording was played of a territorial male’s song. Based on the measurements of distance before and after playback, it looks like there isn’t really a pattern of whether or not sex plays a role in an individual bird’s movement after a territorial playback.
lmmtit <- lmer(dist_m ~ sex + (1|sex/id), data = tit)
summary(lmmtit)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dist_m ~ sex + (1 | sex/id)
## Data: tit
##
## REML criterion at convergence: 220.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3793 -0.5409 -0.0890 0.4035 1.6695
##
## Random effects:
## Groups Name Variance Std.Dev.
## id:sex (Intercept) 258.264 16.071
## sex (Intercept) 4.187 2.046
## Residual 99.407 9.970
## Number of obs: 28, groups: id:sex, 14; sex, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 28.838 6.533 14.604 4.414 0.000534 ***
## sex2 -2.629 9.909 14.236 -0.265 0.794568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sex2 -0.659
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Hessian is numerically singular: parameters are not uniquely determinedUsing a linear mixed model to see if there’s an interaction between the sex of the bird and distance of detection, Because the design is perfectly balanced, we can obtain the variance component estimates using mean squares from a conventional ANOVA. But this is confusing because df are different than what we might expect (there were only 12 mosquitos in the experiment).
anova(lm(dist_m ~ sex/id, data = tit))## Analysis of Variance Table
##
## Response: dist_m
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 47.4 47.40 0.1309 0.7207
## sex:id 2 90.6 45.31 0.1251 0.8830
## Residuals 24 8692.3 362.18